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ABSTRACT 

In this paper we outline methodology to efficiently simulate (jump) diffusion bridge sample paths without 
discretisation error. We achieve this by considering the simulation of conditioned (jump) diffusion bridge 
sample paths in light of recent work developing a mathematical framework for simulating finite dimensional 
sample path skeletons (which flexibly characterise the entirety of sample paths). 

1 INTRODUCTION 

Diffusions and jump diffusions arc an important class of stochastic processes widely used to model 
phenomena in a broad range of application areas, such as economics and finance (Black and Scholes 
1973) and the life sciences (Golightly and Wilkinson 2006). Diffusions arc also widely used throughout 
computational statistics as their simulation underpins a broad class of highly efficient Markov chain Monte 
Carlo algorithms (Roberts and Tweedie 1996). A jump diffusion V : R—>-R is a Markov process which can 
be defined as the solution to a stochastic differential equation (SDE) of the form (denoting V, := I i m v j 7 V s ), 

dV t = j3(V,.)df + ff(V,.)dW f + d/Vo = v G R, f £ [0, T], (1) 

where /3 : R—>-R and a : R—>-R + denote the (instantaneous) drift and diffusion coefficients respectively, 
W, is a standard Brownian Motion and jf‘‘ denotes a compound Poisson process. is parameterised 
with (finite) jump intensity X : R-»R + and jump size coefficient /./ : RoR with jumps distributed with 
density f\ L . All coefficients arc themselves (typically) dependent on V, and regularity conditions arc assumed 
to hold to ensure the existence of a unique non-explosive weak solution (see (0ksendal and Sulem 2004)). 

We may naturally be interested in simulating sample paths from the measure on the path space induced 
by (1), which we denote by T ( j T . Clearly this is non trivial as sample paths arc infinite dimensional 
random variables (and so at most we can simulate some finite dimensional subset of sample paths) and, 
as a closed form representation of the transition density of (1) will be typically unavailable, we may need 
to resort to time discretisation (Kloeden and Platen 1992) which results in the introduction of error. To 
address these challenges a class of rejection-sampling based algorithms (so called Exact Algorithms as they 
avoid the introduction of error) have been developed to simulate a broad range of diffusions (Beskos and 
Roberts 2005, Beskos, Papaspiliopoulos, and Roberts 2008, Chen and Huang 2013, Jenkins 2013, Jenkins 
and Spano 2014) and jump diffusions (Casella and Roberts 2010, Goncalvcs and Roberts 2013, Pollock, 
Johansen, and Roberts 2015) by means of simulating from an equivalent measure FJ, T . 

In this paper we construct exact algorithms to tackle the related problem of simulating conditioned 
jump diffusion sample paths, which can be represented as the solution to an SDE of the following form, 

dV, = j8 (V f _) d it + a(V t .) d W t + dJ^ , U 0 = v G R, V T = w G R, t G [0, T] . (2) 

A conditioned jump diffusion {or jump diffusion bridge ) is simply a diffusion which in addition to having 
a given start point is also conditioned to have some specified end point. For the purposes of this paper 



we restrict our attention to univariate diffusions and impose a number of additional conditions on the 
coefficients of (1,2) (as detailed in Section 2). 

As in (1), we arc interested in simulating sample paths from the measure induced by (2), denoted 
which (as outlined in (Pollock 2013, Gonqalves and Roberts 2013)) can be achieved by constructing an 
equivalent measure Pq ^ from which sample paths can be drawn. There are two key complications when 
constructing an exact algorithm to simulate conditioned jump diffusions which are not present in simulating 
unconditioned jump diffusions (Pollock 2013). Firstly, construction of an appropriate equivalent measure 
Pq w t is more difficult. Secondly, the computational cost of simulating conditioned (jump) diffusions does 
not necessarily scale linearly as a function of the time interval in which it has to be simulated over, and 
so exact algorithms can be rendered computationally infeasible for particular applications. 

In this paper we outline methodology to simulate conditioned jump diffusion sample paths, employing 
strategies to accelerate acceptance and rejection of proposal sample paths and reduce overall computational 
cost. We achieve this by considering the simulation of conditioned jump diffusions in light of recent work 
developing a mathematical framework for simulating diffusion sample path skeletons (characterising the 
entirety of sample paths), and the extension of exact algorithms to Adaptive Exact Algorithms (which enable 
the simulation of lower dimensional skeletons) (Pollock, Johansen, and Roberts 2015). 

This paper is organised as follows: In Section 2 we introduce the key concepts, framework and conditions 
imposed in establishing the results presented in this paper. In Section 3 we introduce more formally exact 
algorithms and introduce a novel adaptive exact algorithm for simulating conditioned diffusions. Finally, 
in Section 4 we extend our approach to simulating conditioned jump diffusions. 

2 PRELIMINARIES 

In (Pollock, Johansen, and Roberts 2015) a framework for constructing exact algorithms was established in 
which entire (jump) diffusion sample paths could be represented by means of simulating a finite dimensional 
skeleton , guided by three key principles. The skeleton typically comprises a layer constraining the sample 
path. In this section we will begin by reviewing these definitions and principles for exact algorithms below, 
and then outline the notation and conditions imposed to establish the results in this paper. 

Definition 1 (Skeleton) A skeleton (<5A) is a finite dimensional representation of a diffusion sample path 
( [V ~ Ty'jp), that can be simulated without any approximation error by means of a proposal sample path 

d T v ’ w 

drawn from an equivalent proposal measure (P,jT) and accepted with probability proportional to (F), 

' cUt 0,7’ 

which is sufficient to restore the sample path at any finite collection of time points exactly with finite 
computation where V\SA ~ 5?. 

Definition 2 (Layer) A layer R(V), is a function of a diffusion sample path V ~ Pq^ which determines 
the compact interval to which any particular sample path V(oj) is constrained. 

Principle 1 (Layer Construction) The path space of the process of interest, can be partitioned and the layer 
to which a proposal sample path belongs can be unbiasedly simulated, R(V) ~ 2% := P ( 'j^. oR~ l . 

Principle 2 (Proposal Exactness) Conditional on Vq, Vt and R(V), we can simulate any finite collection 
of intermediate points of the trajectory of the proposal diffusion exactly, V ~ P ( j‘ l | R \ lKlV]] . 

Principle 3 (Path Restoration) Any finite collection of intermediate (inference) points, conditional on the 
skeleton, can be simulated exactly, V tl ,..., V tn ~ P ( jj. SA. 

To present our work in some generality we assume Conditions 1-6 hold. A fuller discussion of the 
conditions imposed can be found in (Pollock 2013, §1.3, §4.2 & §5.4). 

Condition 1 (Solutions) The coefficients of (1,2) are sufficiently regular to ensure the existence of a unique, 
non-explosive, weak solution. 



2 PRELIMINARIES 


Condition 2 (Continuity) The drift coefficient /3 G C 1 . The volatility coefficient a £ C 2 and is strictly 
positive. 

Condition 3 (Growth Bound) We have that 3K > 0 such that \P(x)\ 2 + |cr( jc)| 2 < K(\ + \x\ 2 ) \/x € R. 
Condition 4 (Jump Rate) A is non-negative and there exists a constant A < °o such that A < A. 

Conditions 2 and 3 are sufficient to allow us to transform our SDEs in (1,2) into one with unit volatility 
(letting \j/[ ,..., i/w r denote the jump times in the interval [0,7’], y/o :=0 and y/y, _ ] — := \^n t +\ '■= T, and 
N t :=£,->i 1 {I/// < t} a Poisson jump counting process). As noted in (Art-Sahalia 2008), this transformation 
is typically possible for univariate diffusions and for many multivariate diffusions. 

Result 1 (Lamperti Transform (Kloeden and Platen 1992, Chap. 4.4)) Let r] (V,) =: X, be a transformed 
process, where p (V t ) := 0 I /a(u)du (where v* is an arbitrary element in the state space of V), then by 
applying Ito’s formula for jump diffusions to find dA, we have (where /./, M-,V t .)=M-,r ] -\X t .))), 


dX, = 




We denote the measure induced by the transformed unconditioned jump diffusion (3) as Q,j T (with 
left hand point Xq := x = p(v)), and Q ( j'' r as the measure induced by the transformed conditioned jump 
diffusion (constrained to have end point Xj := y = p(w)). We further denote by W ( j T as the measure 
induced by the following driftless jump diffusion with unit volatility, 

dA, = d W, + d0 S , A 0 = x G R, t G [0, T], (4) 


where 0' 5 is a compound Poisson process with constant finite jump intensity A, jump size coefficient 
8 : R—>R and with jumps distributed with density f§. We denote by J :{] r as the trajectory of a compound 
Poisson process over [0.7j and as the measure induced by (4) where we additionally have Xj = y. 

In order to deploy an exact algorithm we need to establish that the Radon-Nikodym derivative of 
Qq ^ with respect to Wq ’ y T exists (Results 2 and 3) and can be bounded on compact sets (Result 4). In 
order to do so we impose on the coefficients of (3,4) the following final conditions (where we denote by 
A(u ) := 0a(y)dy and set <I>(X S ) := cr(A s )/2 + a'(X s )/2), 

Condition 5 (<4>) There exists a constant <f> > — °o such that d> < (j). 

Condition 6 (x) We have that 3 x < °o such that, 

MM -fv (X W ,X v 0 .e-lMW-Mx*-)] 

A-f s (X ¥ -X Vi J) 

First considering the Radon-Nikodym derivative of (Qq t with respect to WJ T we have, 

Result 2 (Unconditioned Radon-Nikodym derivative (0ksendal and Sulem 2004)) Under Conditions 1-3 
and 6, the Radon-Nikodym derivative of (Qq t with respect to WJ T exists and is given by Girsanov’s 
formula as follows, 

d ^°/ (A) = exp |A(A r ) — A(x) — J (j) ( A s _ ) d^ | • exp | ^ [A(A,_) - A] dsj 

»r U(X v ,_)-f v (X ¥ -X ¥ 0 
U A-f d (X ¥ -X Vi .) 



In the particular case where we have a diffusion (where X = A = 0), we have, 


- [m) 4 >} ■ 

Now considering the Radon-Nikodym derivative of with respect to Wjp we further denote 
by Pr{x,y) := G dy|X 0 = x)/dy and w T (x,y) := ¥ Wqt (X t G dy|X 0 =x)/dy as the transition 

densities of (3) and (4) respectively over the interval of length 7 initialised at Xq = x. 

Result 3 (Conditioned Radon-Nikodym derivative (Dachuna-Castelle and Florens-Zmirou 1986)) Following 
directly from Result 2 we have, 


dOpT , y . _ wr{x,y) d(Q qj 
dW x Q y ’ r p r (x.y) ' dWJj 

with transition density of the following form (by taking expectations with respect to W'J’’ r ), 


Pr(x,y) 


w T (x,y) ■ E w j, 


dQo.r 

dWj) 


Throughout this paper we rely on the fact that upon simulating a path space layer (see Definition 2) 
then Vs G [0,7"] <j>(X s ) is bounded, however this follows directly from the following result. 

Result 4 (Local Boundedness) By Condition 2, a and a' arc bounded on compact sets. In particular, 
suppose 3£, v G IR such that V t G [0,7], X t (co) G [£,v\ 3Lx := L(X(m)) G := U {X(a >)) G 1R such 

that V t G [0,7], 0(X f (m)) G [Lx,Ux\- 


3 EXACT SIMULATION OF CONDITIONED DIFFUSIONS 

In this section we outline exact algorithms to simulate sample path skeletons of diffusion bridges (under 
Conditions 1-5 and following the Lamperti transform (Result 1)) which can be represented as the solution 
to the following SDE, 


dX T = a{X r )dt+dW t , X 0 =xG ]R,X r =y G IR.,/ G [0,7], (5) 

We present two separate exact algorithms to simulate conditioned diffusion sample path skeletons - 
the Conditioned Unbounded Exact Algorithm (CUEA) and the Conditioned Adaptive Unbounded Exact 
Algorithm (CAUEA) (which is a Rao-Blackwellisation of the CUEA requiring less simulation of the sample 
path). The methodology developed in this section is a direct extension of that developed for unconditioned 
diffusions in (Pollock, Johansen, and Roberts 2015) (termed the Unbounded Exact Algorithm and Adaptive 
Unbounded Exact Algorithm respectively), but also serves to introduce the key ideas for when we consider 
the non-trivial extension to the simulation of jump diffusion bridge sample path skeletons in Section 4. 

Exact algorithms arc a class of rejection samplers operating on diffusion path space (introduced by 
(Beskos and Roberts 2005)) in which finite dimensional subsets of sample paths are drawn from (recall, 
the measure induced by (5)) by means of simulating finite dimensional subsets of sample paths from an (easy 
to simulate) equivalent measure with bounded Radon-Nikodym derivative . As established in Section 2, Wj * y T 
is such an equivalent measure (Brownian motion measure, from which finite dimensional subsets of sample 
paths can be drawn without error (see (Pollock 2013, §2.8))). Proceeding as in standard rejection sampling, 

if we draw X ~ Wj ’ y T and accept the sample path (7=1) with probability ^(*):=»dSf(*)G[0,l] 

then (X\I = 1) ~ Qq Now, considering the form of the acceptance probability we have, 
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Theorem 1 (Conditioned Exact Algorithm Acceptance Probability I) (Qq ’ ;y T is equivalent to Wjy with 
Radon-Nikodym derivative : 

d ^x T y (x) °< exp | - ^ (j) (X s ) ds j € [0, e- 0T ] , (6) 

and so we have that, 

Pw% (*) = e* T • exp | - J* $ (X s ) djj € [0,1]. (7) 

Proof. LHS of (6) from Results 2 and 3. RHS of (6) from Condition 5. (7) rearranged from (6). □ 

As remarked in Section 1, it isn’t possible to simulate entire diffusion sample paths (they are infinite 
dimensional) and so it isn’t possible to evaluate the integral in (7). However, it was noted in (Beskos 
and Roberts 2005, Beskos, Papaspiliopoulos, and Roberts 2006, Beskos, Papaspiliopoulos, and Roberts 
2008) (and summarised in Algorithm 1) that by first simulating an auxiliary random variable F ~ F, an 
unbiased estimator of (6) can be constructed and evaluated without having to simulate the entire sample 
path (i.e. F informs us as to which pa its of the sample path to simulate (denoted X 1 ™)). The remainder 
of the sample path (denoted X rem := X \ A'" 1 ) can be simulated as required after acceptance (hence the 
asterisk in Algorithm 1 Step 4) conditional on its skeleton (composed of F, Xq = x, Xj = y and A 1 " 1 ). 

Algorithm 1 Exact Algorithm for Conditioned Diffusions. 

1. Simulate F ~ F. 

2. Simulate X fin ~ F. 

3. With probability /\ v m f (A) accept, else reject and return to Step 1. 

4. * Simulate X rem ~ (X> ,F). 

Now we consider how to construct a suitable finite dimensional random variable F ~ F (while ensuring 
we satisfy Principles 1-3). As noted in Section 2, to simulate a sample path skeleton we will typically 
require a path space layer. This is due to the fact that the method employed to construct F requires upper 
and lower bounds for 0(X[ OT ]) which, as a consequence of Result 4, is provided by a path space layer 
(Ux G IR and Ly G R respectively). As such the first step in simulating F is to partition the path space of 
Wq-'j. into disjoint layers and simulate the layer to which our proposal sample path belongs (see Principle 
1, denoting R := R(X) ~ ^ as the simulated layer). As such we have for all test functions H G Tf, 

P W f r (X)-H(X) =E^E w .v,v|^ P w w(X)-H(X) . 

Conditional on the simulated layer we can represent our acceptance probability as follows, 

P W * 0 ? r CA) = e~( £ ^ r •exp|-^ r (<j»(A,)-L x ) d sj =: e~( L ^ T ■ P W ^ R (X), (8) 

noting that, 

Pw%\r(X) € [e~^- L ^ T , l] C (0,1], (9) 

As noted in (Beskos, Papaspiliopoulos, and Roberts 2006) and with the aid of Figure 1, (X) is 

precisely the probability a Poisson process of intensity I on the graph 'X A := {(x.y) G [0,7’] x [<&,°o): y 




contains no points. This process can be simulated using a Poisson thinning argument, by means of simulating 
a Poisson process of intensity 1 on the larger graph <£p := [0.7"] x [<t>, Ux] D % (which is trivial), computing 
<j) (X) at a finite collection of time points and then determining whether or not any of the points lie in 
With reference to (8) and as noted in (Pollock, Johansen, and Roberts 2015) (and in the related later 
equivalent construction of (Dai 2014)), this approach to simulating an event of probability P W ' > (X) can be 
made computationally more efficient be deploying an accelerated rejection strategy, in which the sample 
path is first rejected with probability 1 — e -( L x-®) T ((z [q. 1), the crosshatched region in Figure 1) and then, 
conditional on not having been rejected, acceptance is determined by simulating an additional event of 
probability PWg y T \R (X) (the vertically hatched region in Figure 1 which can be simulated as per P W*/ T ( X )’ 

but with the alternate graphs of ■= { (-M’) € [0.7"] x [Lx,°°) : y < 0 (x) } and c £p := [0. T] x [L x JJx])- 
The critical observation in these approaches is that the acceptance probability can be evaluated using only 
a finite dimensional realisation of the sample path, X fin . The above argument is stated more formally in 
Theorem 2, with Algorithm 2 detailing how to implement this strategy to simulate sample path skeletons. 


Ux 


* 


Lx 



Time 


Figure 1: Example trajectory of <j) (X) where X ~ 




MX)- 


Theorem 2 (Conditioned Exact Algorithm Acceptance Probability II (Pollock, Johansen, and Roberts 

2015, §3.1)) Letting K# be the law of K ~ Poi((LA — Lx)T), HJfC the distribution of (<§i,.<§ K ) ~ U[0, T] 
we have, 


P w ^(X)=e-^~^ T - E Ks 




±( Ux-<KXz l ) \ 

Ux-Lx ) 


X 


X 


The computational cost of the CUEA is intrinsically linked to the area of the graph %, and so we 
naturally want to choose or construct the graph ( Xp to occupy as small an area as possible. It was noted 
in (Pollock, Johansen, and Roberts 2015, §3.2) that Algorithm 2 Step 3a could be equivalently performed 
by means of simulating exponential random variables. We could for instance set qo = 0 and iteratively set 
%i = %i -1 + Q where £,■ Exp (LA — Lx) while < 7", or in any other convenient order provided we 
have coverage of the interval [0,7"]. The key idea in (Pollock, Johansen, and Roberts 2015, §3.2) is to 
use this iterative simulation of the sample path to construct an Adaptive Exact Algorithm in which we find 
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Algorithm 2 Conditioned Unbounded Exact Algorithm (CUEA). 

1. Simulate layer information R as per (Pollock, Johansen, and Roberts 2015, §7.1). 

2. With probability (1 — exp { — (Lx - T>)T \) reject path and return to Step 1. 

3. Simulate skeleton points (Xp ,... ,X? ) R. 

(a) Simulate K ~ Poi((t/x — Fx)T) and skeleton times <§i,..., t, K ~ U[0, T], 


4. 


5. 


W $ 


R as per (Pollock, Johansen, and 


(b) Simulate sample path at skeleton times Xt ,... ,Xr 
Roberts 2015, §7.1). 

With probability flMi [(Ux — <j>(Xp.)') /(Ux — Lx)], accept path, else reject and return to Step 1. 


* Simulate X' 




£ Vw $_ 1i6 


R as per (Pollock, Johansen, and Roberts 2015, §3.1). 


refined upper and lower bounds for segments of 0(A[o.r])’ and hence accelerate the acceptance or rejection 
of the sample path (in essence find a smaller graph 'X P to conduct the remainder of the simulation). This 
approach is well suited to simulating conditioned diffusion sample paths as, as noted in Section 1, over 
long time intervals the computational cost for employing an exact algorithm for conditioned diffusions can 
be infeasible (the bounds on the path space layer are less tight and hence the graph 'X P is larger). 

As discussed in (Pollock, Johansen, and Roberts 2015, §3.2), the most computationally efficient order of 
simulating the exponential random variables is iteratively emanating from the centre of uncovered intervals 
(where there is the opportunity to learn most about the extent to which the sample path oscillates). In 
particular, beginning at the interval mid-point ( T /2), we can find the skeletal point closest to the mid-point 
by simulating T ~ Exp(2[t/x — Lx]) and setting the skeletal point (£,') to be with equal probability either 
772 —Tor 772 +r. Halting our simulation of (9) at this point we arrive at (10) where we have decomposed 
our acceptance probability into the product of three probabilities associated with three disjoint sub-intervals 
(conditional on <7 6 [0, T), we have [0, T] = [0, L/2 — t] l+l [T/2 — T, L/2 + t] l+l [ T /2 + T, T]). If we consider 
the evaluation of each successively we need only continue to the next (and expend computation) conditional 
on the previous being accepted (i.e. we have an accelerated rejection strategy). We begin by evaluating 
the computationally cheap expectation in (10) (which is with respect to u ~ U[0,1]), before proceeding to 
the acceptance probabilities for the left and right sub-intervals, each of which has the same form as (9). 

Ux-MX?) 

Ux-Lx 

/o r/2 " T [0(Xv 

if <^[0,L]. 

Considering in isolation the acceptance probability corresponding to the interval [0, L/2 — t] in (10), we 
can now find new layer information (7?^ ) which more tightly bounds the sample path and so find tighter 

bounds for <7^7^7 (denoted ^ and l|*). As such the acceptance probability can be re-written, 

T T 

exp {~J o L x ]dsj = exp |-(l| ),s ] -L x ) • ( 772 - t )} -exp {-^~ [<$>(X S )-L^ 'jdsj. (11) 

The form of (11) now coincides with (8) and so can be evaluated using the same procedure outlined above. 
Iterating this procedure until the entire sample path is accepted or rejected results in the Conditioned Adaptive 
Unbounded Exact Algorithm (CAUEA) presented in Algorithm 3. In Algorithm 3 we use the following 
notation: n denotes the set comprising information required to evaluate the acceptance probability for each 
interval still to be estimated, n := {n(f)}J^. Each Hff) comprises information regarding the time interval it 
applies to [s(n(i)),f(Il(i))], the sample path at known points at either side of this interval (x{\\ii)) :=X? , 


r W; 


™\rxX x ) ~ j 


E 1 

I 

•exp | 
1 , 


u < 


Xp, 


) Lx] ds J T /2 +T [(j)(X s ) Lx] ds| 


, if Q £ [0, T], 


( 10 ) 




y(II(i)) :=X^^) and the associated layer and induced bounds on (U™ 1> and L^^), noting that 

s < s <t <t. We further denote 2m(H(i)) := [j(n(i)) + f(II(i))], 2d{Yl(i)) := [f(II(/)) — s(n(/))]. 


Algorithm 3 Conditioned Adaptive Unbounded Exact Algorithm (CAUEA). 


1. 

2 . 

3. 

4. 


Simulate layer information Rx ~ M as per (Pollock, Johansen, and Roberts 2015, §8.1), setting 
U:={E}:={{[0,T},X 0 ,X t ,R x }} and k = 0. 

With probability (1 — exp { — (Lx — <b)r}) reject path and return to Step 1. 

Set 2 = 11(1). 

Simulate z ~ Exp (2 [U x — L x ])- If z > d( 2) then set n := n\2 else, 

(a) Set K = fc + 1 and with probability 1 /2 set <§(. = m{ 2) — z else E,' K = m( 2) + z. 

(b) Simulate AT ~ R x as per (Pollock, Johansen, and Roberts 2015, §8.2). 

(c) With probability (1 — [U x — 0 {Xr> K )\/{U x —L x \) rc j ccl sample path and return to Step 1. 


(d) Simulate new layer information R X ~'^ K 
Johansen, and Roberts 2015, §8.3 & §8.4). 


and R 




(e) With probability ^ I — exp j — 
and return to Step 1. 

(f) Set n : = nij {[^(2 
If |2| / 0 return to Step 3. 

Define skeletal points ..., as the order statistics of the set {<§{,..., 


m)&) 


+ L 


[&.J(S)] 


- 2 L x 


5. 

6 . 


conditional on Rf as per (Pollock, 
[d(2) — t] | j reject sample path 


m(2)-T],X/,X^,4 (E) ^ju{[m(S)+T,f(2)],Z^,Xf,Rp (E)1 }\2. 


7. * Simulate X rem 








[ 6 - 1 . 6 ] 


as per (Pollock, Johansen, and Roberts 2015, §8.5). 


Accepted sample path skeletons simulated under both the CUEA and CAUEA arc composed of given 
terminal points, skeletal points and layer information and have a form as shown in (12). Both approaches 
satisfy Principles 1-3 (although, the CUEA requires augmentation with additional layer information as per 
(Pollock, Johansen, and Roberts 2015, §3.1)). In Figures 2(a) and 2(b) we present illustrative examples of 
accepted sample path skeletons under the two approaches. 

^CUEA (X) := {(^Co 1 ,*} , ^CAUEA (X) := { } • (12) 


4 EXACT SIMULATION OF CONDITIONED JUMP DIFFUSIONS 

In this section we extend the methodology of Section 3, outlining how to simulate sample path skeletons 
of conditioned jump diffusions (under Conditions 1-6 and following the Lamperti transform (Result 1)) 
which can be represented as the solution to the following SDE (denoting X, := \ \m^ t X s ), 

dX t = a{X t .)dt + dW, + dJ^ v , X 0 =x <E lR,X r =y <E IR, t € [0,L], (13) 

The approach we take in this section in constructing our exact algorithm is based upon the recent methodology 
developed in (Goncalvcs and Roberts 2013). However, we reformulate the exact algorithm presented in 
(Gonsalves and Roberts 2013) to ensure that upon accepting a sample path skeleton then it is possible to 
simulate the sample path at further finite collections of time points (i.e. it satisfies Principles 1-3) and in 
order to employ accelerated rejection strategies to reduce the computational cost of simulation. 

The rejection sampling construction of Section 3 to simulate sample skeletons from Q ( y{, cannot be 
directly employed in the case of conditioned jump diffusions (13) with as the proposal measure, as 
it is not possible to simulate a compound Poisson process conditioned to hit a specified end point. The key 
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(a) Example skeleton output from the 
CUE A (Algorithm 2), Acuea (A), over¬ 
laid with two possible sample path tra¬ 
jectories consistent with the skeleton. 


(b) Example skeleton output from the 
CAUEA (Algorithm 3), Afc AUEA (A), 
overlaid with two possible sample path 
trajectories consistent with the skeleton. 


(c) Example skeleton output from the 
CAUJEA (Algorithm 4), J^caujeaPO, 
overlaid with two possible sample path 
trajectories consistent with the skeleton. 


Figure 2: Comparison of the CUEA, CAUEA and CAUJEA skeleton output. Flatched regions indicate 
layer information, whereas the asterisks indicate skeletal points. 


contribution of (Goncalves and Roberts 2013) was to note that an alternate equivalent measure (denoted 
<G'o ' y T ) can be constructed to ensure the end point is hit. In particular, if a compound Poisson process is 
simulated first (/[o,r]) then, to ensure the end point is hit {Xj = y), a Brownian bridge conditioned to start at 
Xq :=x' =x and end at X' T := / = y — J T can be used as the continuous component in the proposal sample 
path. Considering the supeiposition of the compound Poisson process sample path and the Brownian bridge 
sample path (X t = J t + A/), then the resulting sample path starts and ends at the desired points (Ap = a and 
Xj = y). More formally (Gq ,y T is the measure induced by the following SDE, 

dX, = dZ, + dJr 5 , A 0 = a <E R, X r =y <E R, f E [0,T], (14) 


where Z ~ Wj’^ (where W'Jj is Brownian bridge measure starting at Zo = x and ending at Zj = y' = y — Jr). 

Proceeding as in Section 3, we require the Radon-Nikodym derivative of with respect to 
Theorem 3 (Radon-Nikodym derivative for conditioned jump diffusions (Goncalves and Roberts 2013, 
Lemma 2) (Pollock 2013, Thm. 5.4.1)) is equivalent to Gqj with Radon-Nikodym derivative : 
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Following our exact algorithm construction of Section 3, if we simply draw A ~ Gyj- and accept the 

i/pv A'.J ’ 

sample path (7=1) with probability P G ^(X) := ^ (A) e [0,1], then we have that (A|7 = 1) ~ Q'y y T . 


























































Considering the form of the acceptance probability (by rearrangement of (15)) we have, 
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As in Section 3, by first simulating a finite dimensional auxiliary random variable F ~ F an unbiased 
estimator of (16) can be constructed and evaluated without having to simulate the entire sample path 
(leaving us with a sample path skeleton). In this instance the first step in constructing F is to follow our 
construction of the proposal measure Gq - ^ in (14), and simulate the process / 0 7 - ~ JX (where JX is the 
law of the compound Poisson process component of Gq y). As such we have for all test functions H G XL, 
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Further denoting by W \ JX as the law induced by simulating (XL ,....XL ) ~ Wj 7 we have, 
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(18) 


Note that our acceptance probability P G x, y ( X ) has been decomposed into three separate acceptance prob¬ 
abilities (all of which need to be accepted). This construction leads naturally to an accelerated rejection 
sampling strategy in which we have a sequence of acceptance probabilities and only proceed to evaluate 

the next conditional on acceptance of the current. Pill. y (X) can be evaluated following the simulation of 

Gq',t 

the compound Poisson process (17), and PYl y (X) can be evaluated once the trajectory of the sample path 
at the jump times is simulated (18). This leaves PLX (X) which has the following form. 


Nt+i ( rVi 

(X) = n ^(M-i) . exp |- [0(X,_) +A(X,_)] dsj . (19) 

Noting that between any two jump times with known end points that no further jumps occur and the sample 
path is a Brownian bridge, then each component of (19) can be considered directly using the methodology 
developed Section 3. In particular - , recalling that <j»(Xis bounded on compact sets, AjXr^^.i) G 
[0, A], denoting Rj := R X [ ¥i , y/, ~ XX as the simulated layer (used to compute Uj := U x \ ¥i _ u¥i ] G 1R and 
Li :=L x \ ¥i _ l:¥i ] G H respectively) then we can compute unbiasedly the required acceptance probability in 
finite computation by means of the following theorem. 



4 EXACT SIMULATION OF CONDITIONED JUMP DIFFUSIONS 


Theorem 4 (Conditioned Jump Exact Algorithm Acceptance Probability) Letting be the law of 

«•(/) ~ Poi((A +Ui-Lj) ■ (Wi-Wi- 1)) and the distribution of (<^i,...,«■(<)) ~ U[y//_i, Wi\ we have, 



Vk(‘) A + Ui-Lj ) 


X 



Simulating a finite dimensional proposal sample path as suggested above leads to the Conditioned Unbounded 
Jump Exact Algorithm (CUJEA) (which for conciseness is omitted and can be found in (Pollock 2013, 
Algorithm 5.4.1)). However, incorporating the ideas of the CAUEA of Section 3 (Algorithm 3), leads 
directly to the Conditioned Adaptive Unbounded Jump Exact Algorithm (CAUJEA) presented in Algorithm 
4, outputting skeletons of the form in (20). In Figure 2(c) we present an illustrative example of an accepted 
CAUJEA sample path skeleton. 
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Algorithm 4 Conditioned Adaptive Unbounded Jump Exact Algorithm (CAUJEA). 

1. Simulate compound Poisson process J'oj] ~ J? as per (Pollock 2013, §2.9.3). 

2. With probability (1 — Pilly {X)) reject path and return to Step 1. 

3. Simulate . Xy ~ Wqj as per (Pollock 2013, §2.8). 

(2) 

4. With probability (1 — Pl~,L y (X)) reject path and return to Step 1. 

5. For i in 1 to ( Nt + 1), 

(a) Simulate initial layer information P, ~ & as per (Pollock, Johansen, and Roberts 2015, §8.1), 

setting n := {2} := and ff ; - = 0. 

(b) With probability (1 — exp{ — (L x ^^ — 4>) • (t//,• — !//,■_i)}) reject path and return to Step 1. 

(c) Set S = 11(1). 

(d) Simulate t ~ Exp (2[A + C g — Lj]). If x > d(E) then set n := n\E else, 

i. Set Ki = K, + 1 and with probability 1 /2 set ^' Kj = mz — T else = m, + T. 


ii. Simulate ~ $ 


Rf as per (Pollock, Johansen, and Roberts 2015, §8.2). 


iii. With prob. (1 — [A+l/f— )—A (Ac/. )]/[A+t/x— ££]) reject path and return to Step 1. 

[s(S),^'] [£,' k . /(E)] 

iv. Simulate new layer information R t ' and R i ' conditional on Rj as per (Pollock, 
Johansen, and Roberts 2015, §8.3 & §8.4). 

wo u K-rt (r / 

v. With probability I 1 — exp j 

and return to Step 1. 


X[Yi-uYi] ,y/,-] 


vi. Setn:=nu{[5s,m E }u{[ws + T,fs],X^.,X^, 
(e) If Ini / 0 return to Step 5c. 


[dz — T ]} ) re j ect P at h 
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(f) Define skeletal points <§,• .i,..., <§,■ tK . as the order statistics of the set {^/j,..., <§/ }. 
6. Accept sample path skeleton. 


7. * Simulate X rem 


Nt +l 
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./- 1 


R 






ACKNOWLEDGMENTS 


MP would like to thank Flavio Goncalves, Adam Johansen and Gareth Roberts for stimulating discussion 
on this paper. This work was supported by the EPSRC [grant numbers EP/P50516X/1 and EP/K014463/1]. 

References 

Ait-Sahalia, Y. 2008. “Closed-form likelihood expansions for multivariate diffusions.”. The Annals of 
Statistics 36:906-937. 

Beskos, A., O. Papaspiliopoulos, and G. Roberts. 2006. “Retrospective Exact Simulation of Diffusion 
Sample Paths with Applications.”. Bernoulli 12:1077-1098. 

Beskos, A., O. Papaspiliopoulos, and G. Roberts. 2008. “A Factorisation of Diffusion Measure and Finite 
Sample Path Constructions.”. Methodology and Computing in Applied Probability 10:85-104. 
Beskos, A., and G. Roberts. 2005. “An Exact Simulation of Diffusions.”. Annals of Applied Probability 15 
(4): 2422-2444. 

Black, F., and M. Scholes. 1973. “The Pricing of Options and Corporate Liabilities.”. Journal of Political 
Economy 81 (3): 637-654. 

Casella, B., and G. Roberts. 2010. “Exact Simulation of Jump-Diffusion Processes with Monte Carlo 
Applications.”. Methodology and Computing in Applied Probability 13 (3): 449—473. 

Chen, N., and Z. Huang. 2013. “Localisation and Exact Simulation of Brownian Motion Driven Stochastic 
Differential Equations.”. Mathematics of Operational Research 38:591-616. 

Dachuna-Castelle, D., and D. Florens-Zmirou. 1986. “Estimation of the coefficients of a diffusion from 
discrete observations.”. Stochastics 19:263-284. 

Dai, H. 2014. “Exact simulation for diffusion bridges: an adaptive approach”. Journal of Applied Proba¬ 
bility 51 (2): 346-358. 

Golightly, A., and D. Wilkinson. 2006. “Bayesian sequential inference for nonlinear multivariate diffusions.”. 
Statistics and Computing 16 (4): 323-338. 

Gonqalves, F., and G. Roberts. 2013. “Exact Simulation Problems for Jump-Diffusions”. Methodology and 
Computing in Applied Probability 15:1-24. 

Jenkins, P. 2013. “Exact simulation of the sample paths of a diffusion with a finite entrance boundary”. 
arXiv preprint arXiv:1311.5777. 

Jenkins, P, and D. Spano. 2014. “Exact simulation of the Wright-Fisher diffusion.”. Technical report, 
CRISM, Department of Statistics, University of Warwick. 

Kloeden, P, and E. Platen. 1992. Numerical Solution of Stochastic Differential Equations.. 4th ed. Springer, 
Berlin. 

0ksendal, B., and A. Sulem. 2004. Applied Stochastic Control of Jump Diffusions.. 2nd ed. Springer, 
Berlin. 

Pollock, M. 2013. Some Monte Carlo Methods for Jump Diffusions. Ph. D. thesis, Department of Statistics, 
University of Warwick. 

Pollock, M., A. Johansen, and G. Roberts. 2015. “On the Exact and £-Strong Simulation of (Jump) 
Diffusions.”. Bernoulli. 

Roberts, G., and R. Tweedie. 1996. “Exponential convergence of Langevin distributions and their discrete 
approximations”. Bernoulli :341—363. 

AUTHOR BIOGRAPHY 

MURRAY POLLOCK is a Postdoctoral Research Fellow in Statistics based at the University of Warwick 
working on the EPSRC programme grant “Intractable Likelihood: New Challenges from Modem Applica¬ 
tions (i-like),” held jointly along with Bristol, Lancaster, and Oxford universities. His research interests lie in 
Monte Carlo methodology (particularly MCMC and SMC). His email address is m.pollock@warwick.ac.uk. 


